! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.



PROGRAM CONVNETCDF2GRIB
USE NETCDF
USE GRIB_API
USE LIB_DATES

IMPLICIT NONE

! Namelist variables 
CHARACTER(LEN=300) :: CFNETCDF ! Input netcdf file 
CHARACTER(LEN=300) :: CFSURF   ! input netcdf file surfclim
CHARACTER(LEN=300) :: CFGRIB   ! input template grib file 
CHARACTER(LEN=300) :: CFGRIB_OUT ! output grib file 
CHARACTER(LEN=100) :: CVAR     ! netcdf variable name to process 
LOGICAL            :: LMISS    ! set values in all field if true 
REAL(KIND=JPRB)    :: ZMISS    ! value to set in all field 
LOGICAL            :: LMINMAX  ! if true set min/max bound 
REAL(KIND=JPRB)    :: ZMIN,ZMAX ! min/max value to set 
REAL(KIND=JPRB)    :: DCOEF     ! Coef. to multiply field 
INTEGER            :: INDL     ! vertival leel to read (-1 for none)
INTEGER            :: INDT     ! time index to process (-1 for all, -2 skipping last, -3 skipping last using all grib templates )
INTEGER            :: DATE0    ! BASE DATE  YYYYMMDD
INTEGER            :: TIME0    ! BASE TIME IN HOURS
INTEGER            :: DELH     ! DELT HOURS BETWEEN STEPS
INTEGER            :: DELHT    ! DELT HOURS BETWEEN STEPS in template 
CHARACTER(LEN=5)   :: FCTYPE   ! type of data 'AN', 'FC24', 'FCTOT'
LOGICAL            :: LMASK    ! IF TRUE USE MAKS (IF CFSURF IS NOT NONE)

NAMELIST/NCONF/CFNETCDF,CFSURF,CFGRIB,CFGRIB_OUT,CVAR,LMISS,ZMISS,LMINMAX,ZMIN,ZMAX,&
              DCOEF,INDL,INDT,DATE0,TIME0,DELH,FCTYPE,DELHT

! local variables 
REAL(KIND=JPRB),ALLOCATABLE :: ZINPUT(:),ZOUTPUT(:)

INTEGER*4 :: NCID,VARID,NDIMS,DIMIDS(5),NPCDF,NLEV1,NT1
INTEGER*4 :: NCID_C,VARID_C
INTEGER :: IFILE,IGRIB,IRET,NPGRIB
INTEGER :: IFILEOUT,IGRIB_OUT
INTEGER*4,ALLOCATABLE :: JPOINTS(:)
INTEGER  :: JP,JT,JTR

INTEGER(KIND=JPIM)  :: ZTIME0,CDATE,CTIME,ZTIME
INTEGER(KIND=JPIM)  :: GDATE,GTIME,GSTEP

! set default values 
 CFNETCDF='none'
 CFSURF='none'
 CFGRIB='none'
 CFGRIB_OUT='none'
 CVAR='none'
 LMISS=.FALSE.
 ZMISS=9999.
 LMINMAX=.FALSE.
 ZMIN=1.0e-20
 ZMAX=1.0e20
 DCOEF=1.
 INDL=-1
 INDT=-1
 DELH=-999
 DELHT=-999

! read from namelist 
OPEN(10,FILE="input.nam",STATUS="OLD")
READ(10,NCONF)
LMASK=.TRUE.
IF (TRIM(CFSURF) .EQ. 'none'  ) LMASK=.FALSE.

PRINT*,'CFNETCDF',TRIM(CFNETCDF)
PRINT*,'CFSURF',TRIM(CFSURF)
PRINT*,'LMASK',LMASK
PRINT*,'CFGRIB',TRIM(CFGRIB)
PRINT*,'CFGRIB_OUT',TRIM(CFGRIB_OUT)
PRINT*,'CVAR',TRIM(CVAR)
PRINT*,'LMISS',LMISS
PRINT*,'ZMISS',ZMISS
PRINT*,'LMINMAX',LMINMAX
PRINT*,'ZMIN',ZMIN
PRINT*,'ZMAX',ZMAX
PRINT*,'DCOEF',DCOEF
PRINT*,'INDL',INDL
PRINT*,'INDT',INDT
PRINT*,'DATE0',DATE0
PRINT*,'TIME0',TIME0
PRINT*,'DELH',DELH
PRINT*,'DELHT',DELHT

ZTIME0 = DATE2HH(DATE0,TIME0)

! OPEN INPUT NETCDF FILE 
CALL NCERROR(NF90_OPEN(CFNETCDF,NF90_NOWRITE,NCID),'OPENING NETCDF FILE')
CALL NCERROR(NF90_INQ_VARID(NCID,CVAR,VARID),'FINDING VAR ID')
CALL NCERROR(NF90_INQUIRE_VARIABLE(NCID=NCID,VARID=VARID,NDIMS=NDIMS,DIMIDS=DIMIDS),'FINDING VAR IDS-INFO')
PRINT*,'VARID,NDIMS,DIMIDS',VARID,NDIMS,DIMIDS(1:NDIMS)
CALL NCERROR(NF90_INQUIRE_DIMENSION(NCID,DIMIDS(1),LEN=NPCDF),'GETTING NPOINTS IN NETCDF')
IF (NDIMS==3) THEN
  CALL NCERROR(NF90_INQUIRE_DIMENSION(NCID,DIMIDS(2),LEN=NLEV1),'GETTING NLEVS IN NETCDF')
  CALL NCERROR(NF90_INQUIRE_DIMENSION(NCID,DIMIDS(3),LEN=NT1),'GETTING NTIME IN NETCDF')
ELSE
  NLEV1=0
  CALL NCERROR(NF90_INQUIRE_DIMENSION(NCID,DIMIDS(2),LEN=NT1),'GETTING NTIME IN NETCDF')
ENDIF 
PRINT*,'NPCDF,NLVE1,NT1',NPCDF,NLEV1,NT1

! load mask netcd file
IF (LMASK ) THEN 
  ALLOCATE(JPOINTS(NPCDF))
  CALL NCERROR(NF90_OPEN(CFSURF,NF90_NOWRITE,NCID_C),'OPENING NETCDF CLIM')
  CALL NCERROR(NF90_INQ_VARID(NCID_C,'x',VARID_C),'FINDING VAR ID-X')
  CALL NCERROR(NF90_GET_VAR(NCID_C,VARID_C,JPOINTS(1:NPCDF),(/1/),(/NPCDF/)),'READING NETCDF X')
  CALL NCERROR(NF90_CLOSE(NCID_C),'CLOSING NETCDF CLIM')
ENDIF


!*OPEN template grib file 
IFILE=5
CALL GRIB_OPEN_FILE(IFILE,CFGRIB,'r')
CALL GRIB_NEW_FROM_FILE(IFILE,IGRIB,IRET)
CALL GRIB_GET(IGRIB,'numberOfPoints',NPGRIB)
PRINT*,'NPGRIB',NPGRIB

!* OPEN OUTPUT GRIB FILE 
IFILEOUT=6
CALL GRIB_OPEN_FILE(IFILEOUT,TRIM(CFGRIB_OUT),'w')

!* MAIN TEMPORAL LOOP 
ALLOCATE(ZINPUT(NPCDF))
ALLOCATE(ZOUTPUT(NPGRIB))
IF ( INDT == -2 .OR. INDT==-3 ) NT1 = NT1 -1 !! skip last step (o_gg file )
DO JT=1,NT1

  ! GET DATA FROM NETCDF 
  JTR = JT 
  IF ( INDT >= 1 ) JTR = INDT
  IF ( NDIMS == 2 ) THEN
    CALL NCERROR(NF90_GET_VAR(NCID,VARID,ZINPUT(1:NPCDF),(/1,JTR/),(/NPCDF,1/)),'READING NETCDF VAR')
  ELSE
    CALL NCERROR(NF90_GET_VAR(NCID,VARID,ZINPUT(1:NPCDF),(/1,INDL,JTR/),(/NPCDF,1,1/)),'READING NETCDF VAR')
  ENDIF
  ! APPLY MULTIPLICATIVE FACTOR 
  ZINPUT(:) = ZINPUT*DCOEF
  IF (LMISS) THEN
    ZOUTPUT(1:NPGRIB)=ZMISS
  ELSE
    CALL GRIB_GET(IGRIB,'values',ZOUTPUT)
  ENDIF

  IF (LMASK) THEN
    DO JP=1,NPCDF
      ZOUTPUT(JPOINTS(JP)) = ZINPUT(JP)
    ENDDO
  ELSE
    DO JP=1,NPGRIB
      ZOUTPUT(JP) = ZINPUT(JP)
    ENDDO
  ENDIF 

  IF (LMISS) THEN
    DO JP=1,NPGRIB
      ZOUTPUT(JP) = MAX(ZMIN,MIN(ZMAX,ZOUTPUT(JP)))
    ENDDO
  ENDIF

  ZTIME = ZTIME0 + (JTR-1) * DELH

  !! clone grib file 
  CALL GRIB_CLONE(IGRIB, IGRIB_OUT)
  CALL GRIB_SET(IGRIB_OUT,'values',ZOUTPUT)
  IF (  INDT==-3 )  THEN
     CALL HH2DATE(ZTIME,.FALSE.,CDATE,CTIME)
     IF ( MOD((JTR-1)*DELH,DELHT) .NE. 0 ) THEN ! we have to update the time info here 
        CALL GRIB_SET(IGRIB_OUT,'dataDate',CDATE)
        CALL GRIB_SET(IGRIB_OUT,'dataTime',CTIME*100)
     ENDIF
  ELSE
    IF ( FCTYPE(1:2) == 'AN') THEN
      CALL HH2DATE(ZTIME,.FALSE.,CDATE,CTIME)
      CALL GRIB_SET(IGRIB_OUT,'dataDate',CDATE)
      CALL GRIB_SET(IGRIB_OUT,'dataTime',CTIME*100)
    ELSEIF ( FCTYPE(1:4) == 'FC24') THEN
      CALL HH2DATE(ZTIME,.TRUE.,CDATE,CTIME)
      CALL GRIB_SET(IGRIB_OUT,'dataDate',CDATE)
      CALL GRIB_SET(IGRIB_OUT,'dataTime',0)
      CALL GRIB_SET(IGRIB_OUT,'step',CTIME)
    ELSEIF ( FCTYPE(1:5) == 'FCTOT') THEN
      CALL HH2DATE(ZTIME,.FALSE.,CDATE,CTIME)
      CALL GRIB_SET(IGRIB_OUT,'dataDate',DATE0)
      CALL GRIB_SET(IGRIB_OUT,'dataTime',TIME0*100)
      CALL GRIB_SET(IGRIB_OUT,'step',(JTR-1)*DELH)
    ELSE
      PRINT*, 'FCTYPE not RECOGNIZED: AN, FC24, FCTOT, IN NAMELIST:',FCTYPE
      CALL ABORT()
    ENDIF
  ENDIF
  CALL GRIB_GET(IGRIB_OUT,'dataDate',GDATE)
  CALL GRIB_GET(IGRIB_OUT,'dataTime',GTIME)
  CALL GRIB_GET(IGRIB_OUT,'step',GSTEP)

  CALL GRIB_WRITE(IGRIB_OUT,IFILEOUT)
  CALL GRIB_RELEASE(IGRIB_OUT)
  
!   IF (JT == 1 .OR. JT == NT1 )THEN
    print*,''
    PRINT*,JTR,CDATE,CTIME*100,(JTR-1)*DELH
    PRINT*,JTR,GDATE,GTIME,GSTEP
!   ENDIF 
  IF ( INDT >= 1 ) THEN
    EXIT
  ENDIF
  
  ! UPDATE THE TEMPLATE HERE !
  IF (  INDT==-3 )  THEN
    IF ( MOD((JTR)*DELH,DELHT) == 0 ) THEN  ! then we update the template field 
      CALL GRIB_RELEASE(IGRIB)
      CALL GRIB_NEW_FROM_FILE(IFILE,IGRIB,IRET)
      PRINT*,JTR,'updated template'
    ENDIF
  ENDIF

ENDDO



! CLOSE STUFF
CALL NCERROR(NF90_CLOSE(NCID),'CLOSING NETCDF FILE')
CALL GRIB_RELEASE(IGRIB)
CALL GRIB_CLOSE_FILE(IFILE)
CALL GRIB_CLOSE_FILE(IFILEOUT)


CONTAINS


SUBROUTINE NCERROR(STATUS,STRING)
USE NETCDF
IMPLICIT NONE

INTEGER,INTENT(IN) :: STATUS
CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: STRING

IF ( STATUS /= 0 ) THEN
  PRINT*, TRIM(NF90_STRERROR(STATUS))
  IF( PRESENT(STRING) ) PRINT*,TRIM(STRING)
  PRINT*,'PROGRAM STOP ! '
  STOP
ENDIF

END SUBROUTINE NCERROR

END PROGRAM CONVNETCDF2GRIB
